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Abstract 

Let {Xi,... ,Xn) be multivariate normal, with mean vector pi and covariance 
matrix S, and Sn = e'^^ + • • • + e^". The Laplace transform £.{0) — Ee”®"®” ex 
f exp{—he(a;)} dec is represented as £{6)I{6), where £{0) is given in closed- 
form and I(6) is the error factor (« 1). We obtain £{9) by replacing ho{x) 
with a second order Taylor expansion around its minimiser x*. An algorithm 
for calculating the asymptotic expansion of x* is presented, and it is shown 
that I{0) —>■ 1 as 6^ —>■ oo. A variety of numerical methods for evaluating I{0) 
are discussed, including Monte Carlo with importance sampling and quasi- 
Monte Carlo. Numerical examples (including Laplace transform inversion for 
the density of Sn) are also given. 

Keywords: Lognormal distribution; asymptotics; saddlepoint approximation; 
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1. Introduction 

The lognormal distribution arises in a wide variety of disciplines such as engineering, 
economics, insurance, and finance, and is often employed in modeling across the sci¬ 
ences [3, 13, 16, 22, 23]. It has a natural multivariate version, namely (e'’*^^, •. •, e^") ^ 
LN (/X, S) when {Xi,..., Xn) ^ N {pi, S). In this paper, we consider sums of lognormal 
random variables, Sn =' J- • • • J- e'^”, where the summands exhibit dependence (S 
is non-diagonal), using the notation that Sn ^ SLN(/x, S). Such sums have many 
challenging properties. In particular, there are no closed-form expressions for the 
density f{x) or Laplace transform £{9) of S'„. 
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Models using sums of dependent lognormals are widely applicable, though they are 
particularly important in telecommunications and finance [15, 16]. Indeed, many of the 
approximations for the Laplace transform of sums of independent lognormals originated 
from the wireless communications community [9]. This reflects the significance of 
the SLN distribution within many models, and also that the Laplace transform is of 
intrinsic interest (engineers frequently work in the Laplace domain). In finance, the 
value of a portfolio (e.g. a collection of stocks) is SLN distributed when using the 
assumptions of the common Black-Scholes framework. Thus the SLN distribution is 
central to the pricing of certain options (e.g., Asian and basket) [27]. Also, financial 
risk managers require estimates of f{x) across x G (0,E[S'„]) to estimate risk measures 
such as value-at-risk or expected shortfall. Estimation of this kind has long been a legal 
requirement for many large banks, due to the Basel series of regulations (particularly, 
Basel II and Basel III), so in this context approximating C{9) is useful as a vehicle 
for computing the density f{x) or the c.d.f. These issues are carefully explained in 
[14], [17], and the new Chapter 1 in the recently revised volume of McNeil et al. [26]. 
Comprehensive surveys of applications and numerical methods for the LN and SLN 
distributions are in [21, 5, 6]. 

There exist many approximations to the density of the SLN distribution. Many 
approximations work from the premise that a sum of lognormals can be accurately 
approximated by a single lognormal [10], that is Sn ^ L where L ~ LN(/xl,(t£). We 
refer to this approach as the SLN « LN approximation. Some well-known SLN « 
LN approximations are the Wilkinson-Fenton [18] and Schwartz-Yeh [28] approaches. 
These were originally specified for sums of independent lognormals, but have since been 
generalised to the dependent case [2]. A more recent procedure (for the independent 
case) is the minimax approximation of Beaulieu and Xie [11] calculating the values of 
fiL and which minimise the maximum difference between the densities of S'„ and L. 
However, [11] concludes that the approach is inaccurate in large dimensions or when 
the Xi have signihcantly different means or standard deviations. Finally, Beaulieu and 
Rajwani [10] describe a family of functions which mimic the characteristics of the SLN 
distribution function (in the independent case) with some success, i.e., high accuracy 
and closed-form expressions. 

Another related avenue of research focuses on the asymptotic behaviour of f{x) 
in the tails. First, Asmussen and Rojas-Nandayapa [7] characterised the right tail 
asymptotics. Next, Gao et al. [19] gave the asymptotic form of the left tail for 
n = 2. Gulisashvili and Tankov [21] then provided the left tail asymptotics for linear 
combinations of n > 2 lognormal variables. Yet these asymptotic forms cannot be used 
to approximate f(x) with precision; to quote [21, p. 29], “these formulas are not valid 
for X > 1 and in practice have very poor accuracy unless x is much smaller than one”. 
Similar numerical experience is reported in Asmussen et al. [6]. 

The approach taken here is via the Laplace transform. Accurate estimates for the 
Laplace transform can be numerically inverted to supply accurate density estimates. 
Asmussen et al. [5, 6] outline a framework to estimate C{0) for n = 1 using a modified 
saddlepoint approximation. In their work, the transform is decomposed into C{0) = 
C{9)I{9), where C{9) has an explicit form and an efficient Monte Garlo estimator is 
given for I(9). 
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This paper generalises the approach of [5, 6] to arbitrary n and dependence. The 
defining integral for the Laplace transform of Sn is 


m 


A/(27T)”det(S) Jb" 


exp 


n 


2=1 


- Dx > da; 

2 i 


( 1 ) 


where D '^= Yi ^ (assuming S to be positive definite so D is well-defined). Write 
the integrand as exp{ —ft,e(a:)}. The idea is then to provide an approximation L(9) by 
replacing hg{x) by a second order Taylor expansion around its minimiser x*. Whereas 
the minimiser x* has a simple expression in terms of the Lambert W function when 
n = 1, as in [5, 6], the situation is much more complex when n > 1. As one of our main 
results we give a limit result for x* as 0 —>■ oo. Further, it is shown that the remainder 
I{9) in the representation C{9) = C{9)I{9) goes to 1, a discussion of efficient Monte 
Carlo estimators of I{9) follows, and numerical results showing the errors of our C{9) 
and (numerically inverted) f{x) estimators are given. The paper concludes with an 
informal discussion regarding estimation of the SLN distribution function F{x), and 
some closing remarks. 


2. Approximating the Laplace transform 

Although the definition (1) makes sense for all 0 G C with 5ft(0) > 0 (we denote this 
set as C+), we will restrict the focus to 0 G (0, oo). Of particular interest are the terms 
in the exponent, which in vector form (see Remark 2.1 below) are 

hg{x) =* 0(e'^)^e‘^ + ^x^Dx . 

An approximation of simple form to C{9) — written as C{9) — is available if hg{x) is 
replaced by a second order Taylor expansion. The expansion is given in the proposition 
below. 


Remark 2.1. On vector notation. All vectors are considered column vectors. Func¬ 
tions applied elementwise to vectors are written in boldface, such as e® A (e®i,..., e^")^ 
and log a; “A (loga;i,... ,logXn)^. If a vector is to be elementwise raised to a common 
power, then the power will be boldface, as in “A (^x^, ..., x^)^. The notation x oy 
denotes elementwise multiplication of vectors. The function diag(-) converts vectors to 
matrices and vice versa, like the MATLAB function. o 

Proposition 2.1. The second order Taylor expansion of hg{x) about its unique min¬ 
imiser X* is 

— ^1 — Dx* -I- -(x — x*)~''(A -I- D)(x — X*) 
where A Af 0diag(e'^+“^ ). 

Proof. As hg{x) is strictly convex, a unique minimum exists. Since '\/h${x*) = 0, 
the linear term vanishes in the Taylor expansion, so we have 

hg{x) « hg{x*) ~\- i(x — x*)~^H{x — X*) 
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where H is defined as the Hessian {d‘^hg{x)/dxi dxj) evaluated at x*. To find the 
value of H, we just take derivatives: 

\/hff{x) = de>^+^+Dx, iT = A + D = £)(/ +SA). 

Since A and D are both positive definite, so is H. Also, Vhg{x*) = 0 gives 

—= Dx* which implies —= 1^Dx*. (2) 

Therefore the expansion becomes 


hg{x) 


-l^Dx* + ^{x*yDx* + ^{x- ®*)^(A + D){x - X*) 
— ^1 — 2 ®*) Dx* + -(x — a;*)^(A + D){x — x*). 


□ 


This expansion allows C{0) to be approximated as a constant factor exp{—/i 0 (a;*)} 
times the integral over a normal density (with inverse covariance A + £>), which leads 
to 

C{e)^C(e)= , ^ =exp|fl--a;*')^£)a;4 . 

y/det(/ + SA) \V 2 / j 

We need a suitable error or correction term in order to assess the accuracy of this 
approximation, so we will decompose the original integral (1) into C{9) = C{9)I{9). In 
the integral of (1) change variables such that x = x* + + SA)”^/^y. Then by 

applying ( 2 ), multiplying by expjl^iDa;* — 1^Dx*}, and rearranging, we arrive at 


m 


^(27T)"det(J + SA) Jb- '■ 

- + SA)-i/ 2 y)Tr>(Si/ 2 (J + SA)-i/ 2 y)} dy 

mm 


where 

jRrr a/(27T)" ^ V 

- + SA)-i/2y) - iyT(7 + SA)-iy} dy. 

This form may not be particularly elegant. However, it can be rewritten in ways more 
convenient for Monte Carlo estimation. 


Proposition 2.2. We have that 

I{9)=E + 'EA)-^/‘^Z) = ^/det{I + SA) E 

where 

g(u) exp I (a;*)^£)(e“ — 1 — m) + itt^SA(J + SA)~^m 
v{u) =‘^ exp {(a7*)^I>(e“ — 1 — m)} , 


( 4 ) 


and Z ~ N (0, J). 
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Proof. To show that I (9) can be written as the first expectation in (4), replace 
in (3) with 

- - ^y^-EAil + SA)~iy + iy^SA(I + SA)-iy 

= -\y^iy + + SA)-iy. 

To prove I(9) equals the second expectation of (4), change variables in (3) to 2 : = 

(J+SA)-l/2y, SO 

I{9) = A/det(J + SA) f ^ e^p{ix*)^D(e^"^^^-l-'E^/‘^z)-lz^Iz}dz. 

Jr’' v (27t)” 1 / 2 J 

(5) 

□ 

Remark 2.2. When n = 1, S = cr^, and y = 0, (5) becomes 

1(e) = VTTe^ “p {-1 - »-) - ■ 

This can be simplified using the Lambert IT function, denoted W(-), which is defined 
as the solution to the equation = 2 [12]. With this we have x* = —yV(9a^). 

Also, we can manipulate 

\/l + 9a‘^e^* = \/l — X* = 1/1 + >V(0cr^) 

so I{9) becomes 

I{9) = ^l + W{9a^) £ ^exp{-^^^(e-^-l-a2)-iz2|d^ 
which coincides with the original result of [5] equation (2.3). o 

3. Asymptotic behaviour of the minimiser x* 

We first introduce some notation. For a matrix X, we write Xi^. and X.^i for the fth 
row and column. Denote the row sums of ZD as a = (oi,..., a„)^, that is, Ui = Di , 1. 
For sets of indices Di and D 2 , then Xq^^^q^ denotes the submatrix of X containing 
row/column pairs in {(u,u) : u G Di,u G D 2 }. A shorthand is used for iterated 
logarithms: log^ 0 = log0 and log„ 0 = loglog„_;i0 for n > 2 (note that log^, 0 is 
undefined for small or negative 9, however this is no problem as we are considering the 
case 9 —>■ 00 ). 

The approach taken to find x* = (a;!,..., is to set the gradient of hg{x) to 0, 
that is, to solve 

6)e'^+»>* + Dx* = 0 . (6) 

We will show that the asymptotics of the x* are of the form 

n 

X* = X!logj 9- yi Pa + ri{9) 


(7) 
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for some /3 = G c = (ci,...,c„)^ G R", and r(6») = (ri(6»),..., r„(6»))T 

where each ri{9) = o(l). Before giving the general result, we consider the special case 
where all Ui > 0 since this result and its proof are much simpler. 

Proposition 3.1. If all row sums D are positive then the minimiser x* takes the form 
X* = - log 0 + log 2 9 - p.i + log Oi + n (9) (8) 

where ri{9) = (!l(log 2 0/log 0) = o(l) for l<i <n, as 9 ^ oo. 

Proof. Inserting (8) in (6) we find 

+ Dx* = (a log 9) o — a log 9 + a log 2 9 — Dfj, + D log a + Dr (9) = 0 . 

Looking at these equations we see that we must have 

lim sup max (0) = lim inf min (0) = 0 , 
g i 0 i 

and to remove the log 2 1? term the main term of ri{6) has to be — Iog 2 0/log0. This 
gives the result of the proposition. □ 

In the general case where some < 0, the asymptotic form of a:* is different from 
(8) and its derivation is much more intricate. 

Theorem 3.1. There exists a partition of {1,... ,n} into Jq and JL such that for 
i G P+ 

X* = - log 6» + logfc. 9- fii +a+ o{l) 

for some 1 < ki < n. All x* in JL follow the general form of (7). In more detail, 
there exists a partition of iF_ into Jv(l) and JL \ .7L(1), sueh that if i G .7L(1) then 
jdi^i < — 1 and if i € iF-\ .7v(l) then 

Pi,! = —1, Pi,2 = . . . = Pi,ki-1 = 0, Pi,ki < 0 

for some 1 < ki < n. Finally we have, writing subscripts + and — for Jq and JL, that 
x_ = + o(l) where C = —Dz\D ^. The sets Jq, F_, .7^(1) and the eonstants 

Pij, Ci, ki are determined by Algorithm 3.1 below. 

See Remark 3.1 for some further remarks on the role of the signs of the row sums. 


Algorithm 3.1: 

1. Let /3._i be the value of w that minimises Dw over the set {la : Wj < — 1}. It 
will be proved in the appendix that the solution has Di .fi. i < 0 when jdi i = —1 
and Di.f3.^i = 0 when Pi^i < —1. Accordingly, we can partition {1, .. . ,n} into 
the disjoint sets 

j;(l)=0, P,{l) = {i:Dj,.(3.,!<0}, 

J^o(l) = {i : p,,i = -1, A..A.1 = 0}, .F_(l) = {i : A,i < -1}. 
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2. For k = 2,...,n recursively calculate /3.^k as the value of w that minimises 

Dw whilst satisfying 

Wi = 0 for i G J^+{k — 1), Wi = 1 for i G F,{k — 1), 

Wi < 0 for i G T'o{k — 1), Di.w = 0 for i G — 1). 

It will be proved in the appendix that the solution has < 0 for z G 

To{k — 1), = 0 when ffi^k < 0 for z G Fo{k — 1), and at least one element 

of Fo{k — 1) has Di^.ff.^k < 0. This allows us to create a new partition by 

j-^{k) = J-^{k -1) U T,{k- 1 ), 

iF,{k) = {z G To{k- 1) : /3j,fc = 0,Di^.p.^k < 0}, 

^o{k) = {z G Jiik - 1) : /3i,fc = 0,Di^,f3,^k = 0}, 

= J^_{k - 1) u {z G J=o{k - 1) : A.fe < 0}. 

Terminate the loop early ifJli(fc — 1) = 0. 

3. Say = J^+{k) and iF_ = F_{k). For each z G let ti to be the index of the first 
element of Di .l3 which is negative, and we have cz = \og{—Di^.(3.^(J. Determine 
the remaining elements (using the same subscript shorthand introduced above) 

by 

c_ = -Df} D +{c+ - /x+) + . (9) 

Proof of Theorem 3.1. We propose a solution of the form (7) and show that when 
the pij are constructed from Algorithm 3.1, the remainder term r^ is o(l). 

The construction allows us to draw the following conclusions for the x*. Let F+ 
and JL be the sets as defined in Step 3 above. Consider individually the indices which 
terminated in the F+ and in the F_ sets. In the first case, there exists a ki with 
1 < fcz < rz such that 

f-1, J = h 

Pi,j ~ I j ~ ki, 

[O, otherwise, 

Insertion in (6) gives 

^ ^ A.ilogjT»- /x + c + r(6>) I , 

\j=ki-l ) 


and Di^.f3.j = 


0, 1 < j < ki - 1, 

<0, j = ki- 1. 


showing that the remainder is o(l). 

In the second case, with i G F+, 

A,i < and Di^.fi.j =0, 1 < j < rz. 
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or there exists 1 < ki < n such that 

[-1, i = i, 

= < 0, < j < ki, and Di . f3.j = 0 for 1 < j < n . 

[<0, j = ki, 

For this case we find 

+ D,,x* = o(l) + A..r(0), 

again showing that the remainder is o(l). Lastly, to show a;_ in terms of x^, consider 
+ Z>_ _|_a:+ + D x = 0. As = o(l), then we can see that 

x_ =—Dz^_D_^^x^ + o{l) = Cxj^ + o{l). □ 

In some cases above, we have been able to write the constant ct as an expression 
involving D and /x. For example, in Proposition 3.1 we have Cj = logUj, and in 
Theorem 3.1 (9) gives the value of Ci for i G A. We can show a similar result in the 
general case for all i G .7^(1), that is, for all i where x* = — Iog0 + log 2 0 —/Xi + Ci + o(l). 


Say F, .7^.(1) and A '= ia the subscripts below, * and ~ refer to these sets. 
Since D is regular, so is D^ ,^. Say that D = D.,. — and denote the 

corresponding row sums by a = (uj, i G F). 

Corollary 3.1. For all i G F, 


X* = - log 6» + log 2 6» - /Xj + log a* + r* (6) 
where ri(9) = o(l) and > 0 as 0 —>■ oo. 

Proof. Let b = —/3.p. We have 

iGFil)U Foil), 

* \>1, iG71(l), 

Split D according to indices in F, and A, then 


D,.b = 


i G J^.(l) = JH, 
iGFoil) U JL(1) = A. 


D~,tb, + = 0 and D„ ,b„ + > 0 . 

The first equation gives = —Dz\D.^,,b,, and this with the second equation shows 

Db, = Dl = a — > 0 , 

thus D has all row sums positive and c, = log (Db,) = log (a). □ 

There are some simple forms of S which fall into the case where all Oi > 0. These 
include the case where all diagonal elements of S are identical, and all non-diagonal 
elements are identical. Note, by positive definiteness of S we must have at least one row 
sum positive. Also, if Xi,... is an AR(1) process, then the resulting covariance 
matrix would have all > 0. Meanwhile, cases where 3 < 0 are not difficult to 

find. For the case n = 2 with variances cr^ < cr| and correlation p, a simple calculation 
gives that both row sums are positive when p < a-ijai, and one is negative when 
p > u\lai (see Gao et al. [19] for the expansion of /(x) as x I 0 for these cases). We 
now list a couple of examples of asymptotic forms of x* for specific fi and S which 
have non-positive row sums. 



Approximating the Laplace transform of the sum of dependent lognormals 


9 


Example 3.1. Consider p = (—10,0,10)^ and 


0.5 

1 

2 \ 


/ 14 

-2 

-2 

1 

3 

4 

, D = 

-2 

1 

0 

2 

4 

10 / 


V -2 

0 

0.5 


Implementing the algorithm gives that 


xl = - \og9 A- log 2 0 + (10 + log2) + o(l), 

X 2 = —2 log 0 + 2 log 2 0 + (20 + 2 log 2 ) + o(l), 
xl = -4 log 0 + 4 log 2 0 + (40 + 4 log 2) + o(l), 


/ -1 1 0 10.69 \ 

(/3|c-/x)= -2 2 0 21.39 \ , D{f3\c-fi) = 

\ -4 4 0 42.77 / 


—2 ^ 
0 0 0 0 

0 0 0 0 


(where unimportant values oi D{/3 \ c — /j,) are replaced by stars). 
Example 3.2. Consider p = (1,2,3)^ and 


/ 0.4545 

0.4545 

0.4545 \ 


( 3 

-0.9 

0.1 \ 

0.4545 

1.7204 

1.8470 

, £> = 1 

-0.9 

2 

- 1.1 

\ 0.4545 

1.8470 

2.9862 / 


( 0.1 

- 1.1 

1 / 


Implementing the algorithm gives that 

x\ = - log 0 + log 2 0 - 1 + log 2.2 + o(l), 

X 2 = — log 0 + logg 0 — 2 + log 0.79 + o(l), 

Xg = — log0 — 0.1 log 2 0 + 1.1 logg 0 — 3 + C3 + o(l) , 

where C 3 = 0.9 — 0.1 log 2.2 + 1.1 log 0.79, and 


f 

1 

0 - 0.2 \ 

- 2.2 

* 

* * \ 

-1 

0 

1 -2.2 , r>(/3|c-p) = 

0 

-0.79 

* * 

V-1 

- 0.1 

1.1 -2.4 J \ 

0 

0 

0 0 / 


Remark 3.1. The importance of the sign of the row sums of D, as illustrated by 
Proposition 3.1, perplexed us for quite some time. However Gulisashvili and Tankov 
[ 21 ] describe an interesting link between the row sums and the minimum variance 
portfolio. They show that the leading asymptotic term of P(S'„ < x) as x 0 depends 
upon 

= min w , where A {w : Wi = 1, Wi > 0} . 

i 

The i in which Wi > 0 indicate which summands in Sn have the ‘least variance’. These 
summands are asymptotically important in the left tail, as they will struggle the most 
to take very small values. Seen from the viewpoint of modern portfolio theory [25] 
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then the solution w is viewed as the optimal portfolio weights to create the minimum 
variance portfolio. When all Oi > 0 then Wi = cLi/ which represents full 

diversification. However when assets become highly correlated (meaning that some 
D row sums are non-positive) then 3wi = 0, i.e., some assets are ignored. Thus the 
asymptotics are qualitatively different when the signs of the row sums change. The 
exact point where an asset’s optimal weight becomes 0 occurs when Oi = 0, and this 
phase change produces a unique and convoluted asymptotic form. As L{9) as 0 —>■ oo 
is related to P(S'ri < x) as x 10 then the behaviour of x* is explained. o 

For applications we will need to find x* for a large number of 9 numerically. The 
results above give a sensible starting point for an iterative solver, such as Newton- 
Raphson. Another option is based on the following formulation. Let A D —diag(£)) 
and write the defining equation as 

-b diag(D) x* = -Ax* . 


For each row i, all x* are now on the left-hand side. Using properties of the Lambert 
W function we see that 


X* = -W 


9ef*' 


■ exp 



X* 


One can use this to perform a componentwise fixed point iteration as an alternative to 
the Newton-Raphson scheme. 


4. Asymptotic behaviour of I(9) 

In order to discuss I{9) as 9 —>■ oo we will consider it in a form different from 
Section 2. Define cr = diag(A-|-D)~^/^ G (0, oo)” and M IS? diag(cr) (A+D) diag(cr) G 
R”""”. In (3), substitute -b = (cr o z), so 


I{9) = 


exp{ — 

— , exp 

A/(27t)"det(M"i) 



— 1 — cr o 2 ; 



dz . 


( 10 ) 


The limit of this integrand is the density of a multivariate normal distribution, which 
when integrated is 1. To see this, consider the following. As 0 —>■ oo we have ( 7 ^ —>■ oo 
or <Ji — >■ Di^i > 0, so taking £ G (2, oo) means 

0 eW+U y = = o(l). ( 11 ) 

Consider the second exponent of (10). For fixed 2 ;, — 1 — aiZi — \crfzf = 0(crf), 

and since af = o(l) by ( 11 ) we have 

0 (eM+=r*)T^e'^o 2 _i_cro 2 _i(cro 2 ) 2 ^ = 0 ( 1 ). ( 12 ) 

Finally, we consider M as 0 —>■ 00 . Say that |Jq| and assume that these are 

the first indices. We can then write that M —>■ M* diag(/„_^, F) where this F 
is the bottom-right submatrix of size (n — 71+) x (n — n+) of the inverted correlation 
matrix implied by S. The M matrices are positive definite for all 9 G (0, 00 ], thus the 
limiting form of the integrand in ( 10 ) is a nondegenerate multivariate normal density. 
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Proposition 4.1. lim I{0) = 1. 

9—^00 

Proof. We use the dominated convergence theorem. By (12) and the paragraph 
which follows that equation, the exponent of the integrand is bounded by a constant 
gi for ||2;|| < 1, say, and that the exponent is below —(/2||^|| otherwise (52 > 0), for 
6 > 6q, say. The latter comes from the positive definiteness of M*, the convergence of 
M to M* , and the convergence of (12). Next, convexity implies that the exponent is 
bounded by (/2||2:|| for ||2;|| > 1. In total we have the bound 

exp - g2||2:||I{||^||>i}|, 

which is an integrable function. Thus the conditions for dominated convergence are 
satisfied and we can safely switch the limit and integral to obtain 7(0) —>■ 1. □ 

5. Estimators of C(6) and I(9) 

The simplest approach is to numerically integrate the original expression in (1). This 
approach is used as a baseline against which the following estimators are compared (the 
approach can, however, be slow or impossible for large n). The next naive approach is 
to estimate the expectation ]E[e“®'^"] by crude Monte Carlo (CMC). This would involve 
simulating random vectors Xi ,..., Xu LN(p, S), with X^ = {Xr^i, ■ ■ ■, Xr^n), and 
computing 

-l K n 

^cmc{9) = ^ ^ exp I - ^ Xr,i^ ■ 

r=l i—1 

However this estimator is not efficient for large 0, and rare-event simulation techniques 
are required. 

Given the decomposition of C{0) = C{9)I{9), then some more accurate estimators 
can be assessed. Simply using C{9) gives a biased estimator (which is fast and determin¬ 
istic) for the transform, however the bias is decreased by estimating I(9) with Monte 
Carlo integration. Proposition 2.2 gives two probabilistic representations of I{9). We 
expect the CMC estimator of the first — E[(/(S^/^(J-|-SA)“^/^Z)] —to exhibit infinite 
variance as 0 —>■ 00 as this has been proven for n = 1 in [5]. Therefore this estimator 
does not seem promising. The second estimator — A/det(T+^A) E[u(S^/^Z)] — 
can be viewed as the first estimator after importance sampling has been applied, so we 
focus upon this. Taking Zi ,..., Zr N(0, S), then 

Ft 

Fis{9) = iexp Dx*\^eyrp {{x*Y - 1 - Z^)) . 

Many variance reduction techniques can be applied to increase the efficiency of these 
estimators. The effect of including control variates into £is{9) was considered, using 
the control variate (x*)^DZf (note the elementwise square). The variance reduction 
achieved was small considering the large overhead of computing the variates (and their 
expectations) so these results have been omitted. Lastly, we considered an estimator 
based on the Gumbel distribution. Say that G = (Gi,...,G„) is a vector of i.i.d. 
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standard Gumbel random variables, that is, P{Gr < x) = exp{—e“^} for a; G R. Then 
C{0) can be rewritten as an integral over the density of a vector of standard Gumbel 
random variables. This estimator was quite accurate, though it had higher relative 
error and variance than the estimators based on £is(6*) so it too has been excluded 
from the results. 

The final two variance reduction technig^ues investigated were common random 
numbers and quasi-Monte Carlo applied to Cis{9); for a detailed explanation of these 
techniques see [20] or [4]. Both individually achieved significant variance reduction, 
and together provided the best estimator. Specifically, 

r=l 

where Qr '= using $~^(-) as the (elementwise) standard normal inverse 

c.d.f., and where {ui, U 2 ,... } is the n dimensional Sobol sequence started at the same 
point for every 9. Therefore, Cq{9) is deterministic (for a fixed R and 9), and using 
this scheme is therefore a kind of numerical quadrature. More sophisticated adaptive 
quadrature methods could possibly be applied. 

6. Numerical Results 

Relative errors are given for the main estimators of C{9) in the table below. In all 
estimators the smoothing technique of using common random variables is employed, 
and all estimators are compared against numerical integration of the relevant integrals 
to 15 significant digits. 


9 

100 

2,500 

5,000 

7,500 

10,000 

C 

-9.89e-3 

-1.27e-2 

-1.28e-2 

-1.27e-2 

-1.27e-2 


1.29e-2 

* 

* 

* 

* 


3.36e-4 

2.96e-4 

2.57e-4 

2.31e-4 

2.11e-4 

^Q. 

-3.19e-6 

-5.03e-6 

-5.31e-6 

-5.56e-6 

-5.98e-6 


Table 1: Relative error for various approximations of C{0) for = 0, S = [1,0.5; 0.5, 1]. 
The number of Monte Carlo replications R used is 10®. Note: a * indicates that the CMC 
estimator simply gave an estimate of 0. 

Also, the p.d.f. of Sn can be estimated by numerical inversion of the Laplace 
transform. As the approximations of C{9) above are only valid for 9 G (0,oo), not 
9 G C+, then this restricts the options for Laplace transform inversion algorithms. 
The Gaver-Stehfest algorithm [29] and so-called power algorithms [8] can be used. We 
report on the results of using the Gaver-Stehfest algorithm as implemented by Mallet 
[24]. 

Other options for estimating f{x) include numerically integrating the convolution 
equation (typically this is only viable for small n), the conditional Monte Carlo method 
(as in Example 4.3 on page 146 of [4]), and kernel density estimation. The following 
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estimators are reported: the conditional Monte Carlo estimator /cond, f = C ^(£(-)), 
/is = and /q £-\£q(-)). 


X 

0.01 

1 

1.5 

2 

3 

^Cond 

-1.17e-l 

2.20e-2 

3.72e-3 

5.21e-3 

-4.60e-3 

I 

-7.03e-3 

2.56e-2 

1.79e-2 

6.00e-2 

3.82e-2 

/is 

1.94e-3 

1.43e-2 

-6.13e-3 

4.00e-2 

3.68e-3 

Iq 

2.90e-4 

l.lle-2 

-9.04e-3 

3.70e-2 

2.44e-3 


Table 2: Relative errors for estimators of f(x) for = 0 and S = [1, 0.5; 0.5,1]. The number 
of Monte Carlo repetitions for each x is H — 10"^ for /cond, /is and /q. 


The numerically inverted Laplace transforms are surprisingly accurate. Using com¬ 
mon random numbers for the £(0) estimators was necessary, otherwise the inversion 
algorithms became confused by the non-smooth input. The precision of the inversion 
algorithms cannot be arbitrarily increased when using standard double floating point 
arithmetic [1], so the software suite Mathematica was used. Yet this did not solve 
the problem of the Gaver-Stehfest algorithm becoming unstable (and very slow) when 
trying to increase the desired precision. Also, the inversion results became markedly 
poorer when f(x) exhibited high kurtosis (i.e., when det(S) became small). 

7. Closing Remarks 

The estimators above give an accurate, relatively simple, and computationally swift 
method of computing the Laplace transform of the sum of dependent lognormals. We 
have shown that the approximation’s error diminishes to zero (1(0) —t 1) as 0 —>■ oo, 
and that it is still accurate for small values of 0. One can find x* — for each 0 examined 
— using a Newton-Raphson scheme and Section 3 gives an accurate starting value for 
the iterations. 

Appendix A. Remaining steps in the proof of Theorem 3.1 

Proof. First we note that all the minimisations are convex problems and therefore 
have unique solutions. 

For the initial step of the algorithm let w be the solution of the minimisation problem 
and let be the vector with 1 at coordinate i and zero at the other coordinates. Then 
gi(s) = (w + esi )T£) (w + esi) is minimised at e = 0. When lUi < — 1 the vector 
w + esi is in the search set for all e: small. We therefore have g( (0) = 0 which gives 
A ,w = 0. When Wi = —1 the vector w + esi is in the search set only for nonpositive 
values of e. This implies g'i(0) < 0 giving Di .w < 0. 

For the general recursive step we let u = and express in terms of 

u from the equations Di .w = 0, i G iF_ (k—1). The derivative of Dw with respect to 
Ui (i being the index inherited from w) is then 2Di^.w+2(dw / dui)D = 

2Di .w. As above we find that the derivatives of Dw with respect to Ui at the 
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minimising point is zero when Ui < 0 and less than or equal to zero when Ui = 0. 

What is left to prove is that Ji{k) always has at least one element with < 

0. To this end define di = — /3..i and dk = dk-i — for k > 1. From the properties 
of f3 we find 

dj-+{k),k = 0 ; dj^^(k),k = 1 and > 0 ; 

dj^o(k),k = 1 and Dj^^k^dk = 0 ; Dj^_(k)dk = 0 . 

Assume now that Di^,(3,^k+i = 0 for all i € Fo{k). We show that this leads to a 
contradiction. Using the assumption /3.,fc+i has the properties 

I^F+(k),k+i = 0; /3j;(fc),fc+i = 1; 

/3jl)(fc),fc+i < 0 and = 0 ; Dj:_(^k)l^>,k+i = 0 . 

Combining the two displays we have 

Dj^(k)dk = and Djr_(^k)dk = Djr_i^k)l3-,k+i- 

Since dk and /3.,fc+i are identical on Jq(fc — 1) and — 1) the equations reduce to 

( d^^k),k ^ = Do \ ^ ^hereDo = ( \ . 

V dj7_(^k),k J V /3j=L(fc).fc / V Djr_,^k),Fo{k) Djr^l^k),F-(k) J 

Since the matrix Dq is positive definite and since dj^(k),k ^ (^j^(k),k we have reached 
a contradiction. □ 
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